Tarea 4: Bayesian Inference

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

Bienvenides a la cuarta tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

#Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tidyr)
library(purrr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
# Análisis bayesiano
library("rethinking")
## Loading required package: rstan
## Warning: package 'rstan' was built under R version 4.3.3
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 4.3.3
## 
## rstan version 2.32.6 (Stan version 2.32.2)
## 
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## 
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## 
## Attaching package: 'rstan'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## Loading required package: parallel
## Loading required package: dagitty
## Warning: package 'dagitty' was built under R version 4.3.3
## rethinking (Version 2.01)
## 
## Attaching package: 'rethinking'
## 
## The following object is masked from 'package:purrr':
## 
##     map
## 
## The following object is masked from 'package:stats':
## 
##     rstudent

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Pregunta 1:

Un conjunto de carteros aburridos de las mordidas de perros ha decidido realizar un catastro de mordidas recibidas por los empleados de su empresa en un periodo de 8 meses, planeando en base a estos datos realizar inferencia bayesiana. Los datos de las mordidas estas datos por el dataset no+mordidas.csv, en donde cada fila representa las mordidas recibidas por diferentes carteros y las columnas señalan si fue mordido o no el cartero en los meses de estudio (si fue mordido es señalado con un 1, de lo contrario es señalado con un 0). Cabe señalar que un cartero no puede ser mordido mas de una vez al mes, ya que el damnificado recibe licencia de todo un mes tras la mordida, reincorporándose el siguiente mes al trabajo.

dataMordidas <- read.csv("no+mordidas.csv", header = T)
cols <- c("bites_month_1", "bites_month_2", "bites_month_3",
          "bites_month_4", "bites_month_5", "bites_month_6",
          "bites_month_7", "bites_month_8")

grid_len <- 100
p_grid <- seq(0, 1, length.out=grid_len)
prior <- dnorm(p_grid, mean = 0.6, sd = 0.05)

posteriors <- c()
counted.months <- c()
p_grids <- rep(p_grid, 4)

for (n in c(1, 2, 4, 8)) {
  subset_data <- dataMordidas[, cols[1:n], drop = F]
  x <- sum(subset_data)
  N <- nrow(subset_data) * n
  
  likelihood <- dbinom(x, size = N, prob = p_grid)
  unstd.posterior <- likelihood * prior
  posterior <- unstd.posterior / sum(unstd.posterior)
  
  posteriors <- c(posteriors, posterior)
  counted.months <- c(counted.months, rep(n, grid_len))
}

post.df <- data.frame("grid"=p_grids, "posterior"=posteriors, "months"=counted.months)

ggplot(post.df) +
  geom_line(aes(x=grid, y=posterior)) +
  facet_grid(months ~ .) +
  labs(title = "Posteriors para diferentes números de meses",
       x = "Proporción de mordidas (p)",
       y = "Posterior")

Respuesta:

A medida que se aumenta la cantidad de meses es posible evidenciar como disminuye la incertidumbre sobre el valor de p dado que se observan colas más pequeñas y una curva más angosta, esto dado que al utilizar más datos la información que proveen estos domina por sobre la prior inicialmente planteada, lo que a su vez demuestra que el cartero planteó que le mordían 6 de cada 10 veces, pero viendo los datos completos se ve que en realidad este valor es más cercano a 4 de 10.

prior <- dunif(p_grid)
posteriors <- c()

for (n in c(1, 2, 4, 8)) {
  subset_data <- dataMordidas[, cols[1:n], drop = F]
  x <- sum(subset_data)
  N <- nrow(subset_data) * n
  
  likelihood <- dbinom(x, size = N, prob = p_grid)
  unstd.posterior <- likelihood * prior
  posterior <- unstd.posterior / sum(unstd.posterior)
  
  posteriors <- c(posteriors, posterior)
}

post.df$posterior2 <- posteriors

ggplot(post.df) +
  geom_line(aes(x=grid, y=posterior, colour="Gaussiana")) +
  geom_line(aes(x=grid, y=posterior2, colour="Uniforme")) +
  scale_color_manual(name="Prior", values=c("Gaussiana"="black", "Uniforme"="red")) +
  facet_grid(months ~ .) +
  labs(title = "Posteriors para diferentes números de meses",
       x = "Proporción de mordidas (p)",
       y = "Posterior")

Respuesta:

Es posible evidenciar que ambas curvas son distintas, con la principal característica siendo que para la prior uniforme, su posterior tan solo refleja la información de los datos, lo que implica una distribución inicial sin preferencias y por ello está más centrada alrededor de 0.4 como se mencionó en el punto anterior, mientras que la gaussiana afecta la forma de la curva haciéndola tender hacia el valor de la media 0.6 que se eligió para esta. Pero a medida que la cantidad de meses aumenta, y por lo tanto, los datos utilizados, la curva con prior gaussiana converge hacia la otra, dado que se disminuyó el impacto que posee su prior inicial al considerar más meses en el análisis.

En base a esto es que se puede decir que la elección inicial de la prior afectará los resultados más notablemente al utilizar pocos datos, pero tras considerar una cantidad mayor esta elección inicial se hará menos relevante, así que al final la importancia de elegir una prior recae en que tan rápido convergerá a la distribución real de los datos.

for (n in c(1, 2, 4, 8)) {
  res <- quap(
    alist(
      W ~ dbinom(L, p), # binomial likelihood
      p ~ dunif(0, 1) # uniform prior
    ),
    data = list(W=sum(dataMordidas[, cols[1:n], drop=F]),L=n*nrow(dataMordidas[, cols[1:n], drop=F]))
  )
  
  df <- precis(res)
  curve(dnorm(x, df$mean, df$sd) / sum(dnorm(x, df$mean, df$sd)), lty=2, add=FALSE,
        main=paste("Posterior usando Laplace (n=", n, ")", sep=""),
        xlab="Proporción de mordidas (p)",
        ylab="Posterior")
}

Respuesta:

Mediante este método se obtienen resultados que difiere principalmente con respecto a la grid-approximation que utilizaba la prior gaussiana, esto debido a que con Laplace se obtuvo una distribución centrada alrededor de 0.4, en contraste con la previa que inicia en 0.6. En cambio, este resultado no se diferencia mucho respecto al caso en que se realiza la grid-approximation con prior uniforme, esto debido a que para el método de Laplace también se utilizó una prior uniforme y los datos poseen un comportamiento simétrico.

¿Cómo puede interpretar los resultados?

grid_len <- 1000
p_grid <- seq(0, 1, length.out=grid_len)
prior <- dunif(p_grid)

W <- sum(dataMordidas)
total <- ncol(dataMordidas) * nrow(dataMordidas)
likelihood <- dbinom(W, size=total, prob=p_grid)

unstd.posterior <- likelihood * prior
posterior <- unstd.posterior / sum(unstd.posterior)

samples <- sample(p_grid, prob=posterior, size=1e4, replace=T)
res <- dens(samples)

Intervalos:

# fracción de valores <0.35
fraction.1 <- sum(posterior[p_grid < 0.35])
fraction.1
## [1] 1.951779e-05
# fracción de valores entre 0.35 y 0.45
fraction.2 <- sum(posterior[p_grid >= 0.35 & p_grid < 0.45])
fraction.2
## [1] 0.9999803
# fracción de valores >0.45
fraction.3 <- sum(posterior[p_grid >= 0.45])
fraction.3
## [1] 2.193479e-07

Respuesta:

En base a los resultados se puede observar como los intervalos (0, 0.35) y (0.45, 1) poseen valores muy cercanos a 0, apoyado visualmente por la densidad de la posterior que se graficó, lo que finalmente indica que una proporción de mordidas en estos intervalos es prácticamente 0, mientras que es casi seguro que la verdadera proporción de mordidas p se encuentre en el intervalo intermedio, es decir, entre 0.35 y 0.45.

Intervalo \(50\%\):

PI(samples, prob=0.5)
##       25%       75% 
## 0.3863864 0.4014014

Intervalo \(75\%\)

PI(samples, prob=0.75)
##       12%       88% 
## 0.3813814 0.4064064

Intervalo \(95\%\)

PI(samples, prob=0.95)
##        3%       98% 
## 0.3723724 0.4154154

Respuesta:

Estos resultados muestran que la proporción verdadera de mordidas está altamente concentrada en torno a valores cercanos a 0.4 dado que el rango es estrecho para todos los casos y además que al aumentar el porcentaje de credibilidad el rango de valores se extiende por muy poco.

Pese a esto, un problema que poseen estos intervalos de credibilidad es que dependen de la distribución de la posterior, y si en particular esta no es simétrica o multimodal entonces estos intervalos pueden omitir regiones significativas de probabilidad fuera del rango seleccionado, aunque para esta experimentación en particular no pareciera ser el caso.

Intervalo \(50\%\):

HPDI(samples, prob=0.5)
##      |0.5      0.5| 
## 0.3863864 0.4004004

Intervalo \(75\%\)

HPDI(samples, prob=0.75)
##     |0.75     0.75| 
## 0.3823824 0.4064064

Intervalo \(95\%\)

HPDI(samples, prob=0.95)
##     |0.95     0.95| 
## 0.3713714 0.4144144

Respuesta:

Se nota que los valores resultantes son muy similares a los obtenidos en la parte anterior, esto indica que la distribución de la posterior es bastante simétrica, esto pues los primeros resultados que se obtuvieron poseen ambas colas con la misma probabilidad asignada, mientras que los HPDI no poseen esta limitación y en cambio solo buscan el rango de valores más estrecho tal que posea en total una cierta probabilidad, de modo que estos dos casos solo son similares cuando la distribución es simétrica.

n_total <- nrow(dataMordidas) * ncol(dataMordidas)
ratio.real <- sum(dataMordidas) / n_total

simulaciones <- rbinom(10000, size=n_total, prob=samples)
ratio.simulados <- simulaciones / n_total

ggplot() +
  geom_density(aes(ratio.simulados)) +
  geom_vline(aes(xintercept=ratio.real), color="red") +
  labs(
    title = "Comparación entre predicciones y proporción real",
    x = "Proporción de mordidas (p)",
    y = "Densidad"
  )

Respuesta:

En base a este gráfico es posible ver como la observación real se haya en el centro de la distribución, mostrando que el modelo se ajusta bien a los datos ya que este es efectivamente un resultado central y probable.

# se eligen el primer y segundo mes para este análisis
mes_1.mordidos <- which(dataMordidas$bites_month_1 == 1)
mordidos_consecutivos.cant <- sum(dataMordidas$bites_month_2[mes_1.mordidos])

prob_mordido2 <- mordidos_consecutivos.cant / length(mes_1.mordidos)
simulaciones <- rbinom(10000, size=length(mes_1.mordidos), prob=prob_mordido2)

ggplot() +
  geom_density(aes(simulaciones)) +
  geom_vline(aes(xintercept = mordidos_consecutivos.cant), color = "red") +
  labs(
    title = "Simulación de carteros mordidos en ambos meses",
    x = "Número de carteros mordidos",
    y = "Densidad"
  )

Respuesta:

Al igual que en el punto anterior, es posible evidenciar que el modelo nuevamente se ajusta de forma correcta a los datos dado que las simulaciones generaron una distribución de datos en la cual el valor observado, en este caso, la cantidad de carteros que fueron mordidos en ambos meses, es central y probable, indicando que describe bien la relación entre las mordidas consecutivas de estos dos meses.

Pregunta 2:

En esta pregunta se trabajará con el dataset “notas.csv” el cual contiene las notas históricas de un curso desconocido. Suponga que los datos vienen de una distribución \(\mathcal{N}(\mu,\sigma^2)\), el objetivo de la pregunta es estudiar el comportamiento de los datos y los posibles valores de \(\mu,\sigma\) mediante técnicas de inferencia Bayesiana.

Usted sabe un dato extra sobre la información, \(\sigma \in [0.5,1.5]\) y, además, tiene una fuerte creencia a que es mas probable encontrar la desviación estándar real entre \([0.5,1]\) que en \((1,1.5]\).

data_notas <- read.csv("notas.csv")
data_notas$Notas <- as.numeric(data_notas$Notas)
# Función para crear likelihood dado mu y sigma
grid_function <- function(mu, sigma) {
  likelihood <- prod(dnorm(data_notas$Notas, mean = mu, sd = sigma))
  return(likelihood)
}

# Valores de la grilla
grid_mu <- seq(1, 7, length.out = 100)
grid_sigma <- seq(0.5, 1.5, length.out = 100)

# Se crea la grilla 2d
data_grid <- expand_grid(grid_mu,grid_sigma)

# Se guarda la likelihod
data_grid$likelihood <- map2(data_grid$grid_mu,data_grid$grid_sigma, grid_function)

# Se transforma el forma de map2 a una columna
data_grid <- unnest(data_grid,cols = c("likelihood"))


# Valores de los priors
prior_mu <- rep(1 / length(grid_mu), length(grid_mu))
prior_sigma <-  ifelse(grid_sigma <= 1,1,1 - (grid_sigma - 1) / 0.5)

# Se crea la grilla 2d de priors
prior <- expand_grid(prior_mu,prior_sigma)

# Se calculan los valores del prior
data_grid$prior <-  map2(prior$prior_mu,prior$prior_sigma, prod)
data_grid <- unnest(data_grid,cols = c("prior"))

# Se calcula el posterior
data_grid$unstd_posterior <-  data_grid$likelihood * data_grid$prior

# Se estandariza el posterior
data_grid$posterior <- data_grid$unstd_posterior/sum(data_grid$unstd_posterior)

# Se ajustan los valores de la posterior para que no sean valores tan pequñeos
data_grid$posterior <- (data_grid$posterior - min(data_grid$posterior))/(max(data_grid$posterior)-min(data_grid$posterior))

Respuesta: Se propuso valores entre 1 y 7 de manera uniforme para el prior de mu, ya que las notas en chile varían entre estos valores, y como no se tiene información del desempeño previo, es razonable pensar que pueden tener una distribución uniforme. Luego, para el prior de sigma se usó la información del enunciado, es decir, está dentro del rango 0.5 a 1.5 con mayor probabilidad entre 0.5 y 1 y menos entre 1 y 1.5, eso determinado en: grid_sigma <= 1,1,1 - (grid_sigma - 1) / 0.5

# Punto de referencia
# Se recomienda cambiar estos valores por unos adecuados que le permitan estudiar
# Los valores de la distribución de mejor manera
valor_x <- mean(data_notas$Notas)
valor_y <- sd(data_notas$Notas)

# Grafico

punto_comparacion <- tibble(x = valor_x, y = valor_y)

plt <- data_grid %>%
  ggplot(aes(x = grid_mu, y = grid_sigma)) +
  geom_raster(aes(fill = posterior),
    interpolate = T
  )+ 
  geom_point(x = valor_x, y = valor_y, size = 1.3,color="white")+
  geom_label(
    data = punto_comparacion, aes(x, y),
    label = "Punto Comparación",
    fill = "green",
    color = "black",
    nudge_y = 0, # Este parametro desplaza la caja por el eje y
    nudge_x = 1 # Este parametro desplaza la caja por el eje x
  )+
  scale_fill_viridis_c() +
  labs(
    title = "Posterior para Mean y Standard Deviation",
    x = expression(mu ["Mean"]),
    y = expression(sigma ["Standar Deviation"])
  ) +
  theme(panel.grid = element_blank())

plt

Respuesta:

El punto de comparación se ubicó en los valores reales de desviación estándar y promedio de las notas. Este se encuentra dentro del área con mayor densidad de la posterior, lo cual indica que el modelo generado concuerda con los datos reales.

Con respecto a la media, su posterior está bastante “concentrada” en un pequeño rango de valores cercano al seis, lo cual evidencia que el modelo estimó de manera precisa el valor (se tuvo la información adecuada). Por otro lado, vemos que los valores de la desviación estándar tienen un poco más de variación que la media (sus probabilidades son mayores a 0 en un mayor rango de valores para el parámetro), pero es común que este parámetro sea un poco más difícil de estimar.

Para un mayor análisis se podrían variar los priors para ver cómo influyen en los resultados y la forma de la posterior.

library(shape)
# Codificamos los datos
x <- 1:length(data_grid$posterior)

# Sampleamos los indices
posterior_samples_aux <- sample(x,size = 1e4, replace = T, prob = data_grid$posterior)

# Obtenemos los verdaderos valores de la sampling distribution
posterior_samples <- data_grid[posterior_samples_aux,]

# Obtenemos solos los valores relevantes para la densidad
df <- data.frame(posterior_samples$grid_mu,posterior_samples$grid_sigma)

# Realizamos las densidades
dens(df)

Respuesta:

Este proceso sirve para aproximar la distribución posterior de los parámetros de un modelo bayesiano utilizando una grilla de los posibles valores para los parámetros. Esto, logrado al hacer sampling de los valores de los parámetros con probabilidad asignada correspondiente. Permite explorar y analizar la distribución posterior de los parámetros, además de evaluar la incertidumbre en las estimaciones de estos.

Del gráfico a la izquierda se puede observar que la densidad varía bastante a lo largo del rango de valores, y la mayor concentración se encuentra alrededor de seis, lo que indica que este es el valor más probable para mu. El derecho tiene mayor densidad cercano a los valores de 0.6 y 0.7, lo que sugiere que los valores de sigma son más probables dentro de ese rango, también se aprecia una incertidumbre relativamente pequeña, ya que en 0.5 y 0.8 la densidad decae significativamente.

Pregunta 3: Regresión Lineal Bayesiana

El objetivo de esta pregunta es introducirlo en la aplicación de una regresión bayesiana. Con esto, buscaremos entender como calcular una regresión bayesiana en base al “motor” aproximación de Laplace, revisando como se comporta la credibilidad de sus predicciones y como la regresión lineal puede llegar a mutar al aplicar una transformación en el vector \(x\). Para responder esta pregunta centre su desarrollo solo en las clases y las funciones que nos brinda la librería rethinking.

Para esta pregunta usaremos el dataset Pokemon.csv, que contiene las características de más de 700 personajes.

poke.df <- read.csv("Pokemon.csv")
poke.df <- na.omit(poke.df[, c("Attack", "Sp..Atk")])
data <- data.frame(Attack = poke.df$Attack, SpAtk = poke.df$Sp..Atk)
data

Respuesta:

library(RcppEigen)
## Warning: package 'RcppEigen' was built under R version 4.3.3
# Modelo Bayesiano
model <- alist(
    SpAtk ~ dnorm(mu, sigma), 
    mu <- a + b * Attack,     
    a ~ dnorm(0, 10),       
    b ~ dnorm(0, 1),
    sigma ~ dunif(0, 50)              
  )

Asumimos que la distribución de SpAtk es normal, al igual que los priors para a y b. Además se establece que la media de la distribución para a y b es igual a 0, y la desviación estándar de a es 10 (mayor espacio a incertidumbre, ya que no se sabe exactamente el promedio de SpAtk y pueden tener un amplio rango de valores) y la de b es 1 (mayor certeza porque se espera que los efectos de Atk sean más limitados). Finalmente, para sigma se establece un prior de una distribución uniforme (asumimos que todos los valores son igualmente probables dentro de un rango) con un límite inferior igual a 0 y superior igual a 50.

model_map <- map(model, 
                 data = data)#map ajusta el modelo usando Laplace approximation
p<-precis(model_map)

a <- coef(model_map)["a"]
b <- coef(model_map)["b"]
sigma <- coef(model_map)["sigma"]

cat(sprintf("Ecuación del modelo:\nSpAtk ~ Normal(mu, sigma)\nmu = %.2f + %.2f * Attack\nsigma = %.2f", 
            a, b, sigma))
## Ecuación del modelo:
## SpAtk ~ Normal(mu, sigma)
## mu = 38.29 + 0.43 * Attack
## sigma = 30.04
Att_seq <- seq(min(data$Attack), max(data$Attack), length.out = 100)
new_data <- list(Attack = Att_seq)

#Predecimos valores de mu link()
mu_pred <- link(model_map, data = new_data)

#Calculamos promedio e intervalos de credibilidad al 95% para la media
mu_mean <- apply(mu_pred, 2, mean)
mu_IC <- apply(mu_pred, 2, PI, prob = 0.95)

#Simulamos valores observados para incluir dispersión residual usando sim()
simulated_values <- sim(model_map, data = new_data)

#Calculamos intervalos de credibilidad al 95% para los valores simulados
sim_IC <- apply(simulated_values, 2, PI, prob = 0.95)
ggplot() +
  geom_point(aes(x = data$Attack, y = data$SpAtk), color = "red") +
  geom_line(aes(x = Att_seq, y = mu_mean), color = "orange") +
  geom_ribbon(aes(x = Att_seq, ymin = mu_IC[1, ], ymax = mu_IC[2, ]), 
              alpha = 0.3, fill = "blue", show.legend = TRUE) +
  geom_ribbon(aes(x = Att_seq, ymin = sim_IC[1, ], ymax = sim_IC[2, ]), 
              alpha = 0.1, fill = "orange", show.legend = TRUE) +
  labs(x = "Attack", y = "Sp. Atk", 
       title = "gráfico de regresión bayesiana")

Lo primero que podemos observar es que la relación entre Attack y Sp.Attack tiende a ser directamente proporcional positiva. Esto sugiere que en promedio, si Attack aumenta, Sp.Atk también. Aquello se aprecia en la distribución de datos reales (puntos rojos) y la regresión generada. Luego, si nos ponemos a analizar el intervalo de credibilidad de 95% para la media (área morada/azul), notamos que tiende a ser bastante estrecha a lo largo de la regresión, lo cual indica que hay una mayor certeza de las predicciones en esta parte. Al comparar con el área para los datos simulados, vemos que esta última es mucho más amplia, lo cual indica una alta variabilidad (hay alta incertidumbre). Con estas observaciones podemos concluir que el modelo es bueno para captar la tendencia promedio entre ambas variables, pero no es el ideal por la alta variabilidad reflejada por el área naranja, ya que indica que Attack no es suficiente para explicar el comportamiento de Sp.Atk.

A work by CC6104